function [lambda, Rdraw, rho] = ...
    BaiWangEstimationSC_SV_new_Extended(schemenumber,y,factors,rho,soft6,Rdrawlast,lambda_priors,trend)
   

n = size(y,2);
nf = size(factors,2);

[R,R1,r] = switch_schemenumber_alt_Extended(schemenumber,nf,n,trend);

stacked_Y =[];
stacked_F =[];
diag_Omega =[];

for i=1:n
              
       if ismember(i,soft6)

          drop = 5; 
           
          Y = y(2+drop+1:end,i);
          X = [];
          for ff = 1:nf
             factor6m = sum([factors(:,ff) mlag2(factors(:,ff),5)],2);
             X = [X factor6m];
           end
           X = X(6:end,:);
       else
       
       drop = 5;
       X=[factors(:,1:nf) mlag2(factors(:,1:nf),5)];
       X = X(5+1:end,:);
       Y = y(2+5+1:end,i);
       end
              
           Y = Y./sqrt(Rdrawlast(i,drop+1:end))';        
           X = X./repmat(sqrt(Rdrawlast(i,drop+1:end))',1,size(X,2));

           Rdraw = eye(n);
       
        Xstar=[];
        Ystar=[];
        Ylag = mlag2(Y,2);
        Ystar = Y +((-rho(i,:))*Ylag')';
        Ystar = Ystar(2+1:end,:);
        
        for jj = 1:size(X,2);
        Xlag = mlag2(X(:,jj),2);
        Xstar(:,jj) = X(:,jj)+((-rho(i,:))*Xlag')';
        end
        Xstar = Xstar(2+1:end,:); 
        
        stacked_Y = [stacked_Y; Ystar];
        stacked_F = blkdiag(stacked_F,Xstar);
        diag_Omega = [diag_Omega; ones(size(Ystar))*Rdraw(i,i)];
        
end


isfinstY = isfinite(stacked_Y);
inv_Omega = spdiags(1./diag_Omega(isfinstY), 0, sum(isfinstY), sum(isfinstY));
stacked_F = sparse(stacked_F(isfinstY,:));
stacked_Y = stacked_Y(isfinstY);


if lambda_priors == 0
    
    product = stacked_F'*inv_Omega*stacked_F;
    lamhatVarhat = InverseOfBlockDiagonalMatrix(product);
    veclamhat = lamhatVarhat*stacked_F'*inv_Omega*stacked_Y;

elseif lambda_priors == 2
    

    matrixOfPriorMeans = [1 zeros(1,5)];

    matrixOfPriorMeans = repmat(reshape(repmat(matrixOfPriorMeans,nf,1),1,[]),n,1);
    veclam_prior = vec(matrixOfPriorMeans');    
    
    tau = 0.2;
    decay = 2;
    
    matrixOfPriorVars = tau./([1:5+1].^decay);   
    matrixOfPriorVars = repmat(reshape(repmat(matrixOfPriorVars,nf,1),1,[]),n,1);
    lamhatVar_prior = diag(vec(matrixOfPriorVars'));    

    lamhatVarhat = pinv(stacked_F'*inv_Omega*stacked_F+pinv(lamhatVar_prior));
    veclamhat =    lamhatVarhat *(stacked_F'*inv_Omega*stacked_Y + (pinv(lamhatVar_prior)*veclam_prior));
end


updatematrix = lamhatVarhat*R'/(R*lamhatVarhat*R');
updatematrix = round(updatematrix*1E14)./1E14;

lambda_r_mean = veclamhat - updatematrix*(R*veclamhat-r);
lambda_r_var = (eye(n*nf*(5+1))-updatematrix*R)*lamhatVarhat*(eye(n*nf*(5+1))-updatematrix*R)';


lambdavec = gmvn(lambda_r_mean,lambda_r_var,logical(R1));

lambda = reshape(lambdavec',nf*(5+1),n)';




